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ABSTRACT 

Reionization is inhomogeneous for two reasons: the dumpiness of the intergalactic medium 
(IGM) and clustering of the discrete ionizing sources. While numerical simulations can in 
principle take both into account, they are at present limited by small box sizes. On the other 
hand, analytic models have only examined the limiting cases of a clumpy IGM (with uniform 
ionizing emissivity) and clustered sources (embedded in a uniform IGM). Here, we present 
the first analytic model that includes both factors. At first, recombinations can be ignored and 
ionized bubbles grow primarily through major mergers, because at any given moment the bub- 
bles have a well-defined characteristic size. As a result, reionization resembles "punctuated 
equilibrium," with a series of well-separated sharp jumps in the ionizing background. These 
features are local effects and do not reflect similar jumps in the global ionized fraction. We 
then combine our bubble model with a simple description of recombinations in the IGM. We 
show that the bubbles grow until recombinations balance ionizations, when their expansion 
abruptly halts. If the IGM density structure is similar to that at moderate redshifts, this lim- 
its the bubble radii to ~ 20 comoving Mpc; however, if the IGM is significantly clumpier 
at higher redshifts (because of minihalo formation, for example), the limit could be much 
smaller. Once a bubble reaches saturation, that region of the universe has for all intents and 
purposes entered the "post-overlap" stage. Because different HII regions saturate over a fi- 
nite time interval, the overlap epoch actually has a finite width. Our model also predicts a 
mean recombination rate several times larger than expected for a uniformly-illuminated IGM. 
This picture naturally explains the substantial large-scale variation in Lyman-series opacity 
along the lines of sight to the known z > 6 quasars. More quasar spectra will shed light on 
the transition between the "bubble-dominated" topology characteristic of reionization and the 
"web-dominated" topology characteristic of the later universe. 
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1 INTRODUCTION 

The reionization of hydrogen in the intergalactic medium (IGM) 
is a landmark event in the history of structure formation: it is the 
epoch at which (radiative) feedback from the first generation of 
sources impacts every baryon in the universe, and thus the first time 
that galaxies can be said to influence the course of structure forma- 
tion throughout the IGM. As such, a good deal of observational 
and theoretical effort has gone into understanding this era. The first 
- and simplest - question to answer is when reionization occurred. 
Currently, the observations paint a complex picture. The Sloan Dig- 
ital Sky Survey (SDSS) has discovered several z > 6 quasars. Each 
shows saturated or nearly saturated absorption in the Ly a forest as 
well as a rapid decrease in the transmission at z ~ 6 iFan et all 
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l200ltlBecker et all200lllwhite et alfeOolllFan et all2O03tl20O4h . 
The appearance of the long-soughT lGTnmj^eTensorr ' ^?65l) ab- 
sorption trough is a crucial constraint on the ionizing background 
of the universe. Although the optical dept h only l imits the global 
neutral fraction to be x m > 10" 3 <Fan et all2002:ISongailal2004 
lOh & Furlanettol2005l) . its r apid evolution may be characteristic o f 
the tail end of reionization <Gnedinll200d: iRazoumov et"ai]f 2002 ) . 
On the other hand, measurements of the polarization of the cos- 
mic microwave background (CMB) by the Wilkinson Microwave 
Anisotropy Probe (WMAP) found a high optical depth to electron 
scattering for CMB photons, which requires reionization to begin 
at z > 15 iKogut et all2 003). albeit with large uncertainties. 

A number of other observations provide additional constraints, 
though their interpretation is even mo re difficult. These in clude: 
the IG M temperature at z < 4 iTheuns et alJl2002l : lHui & Haimai] 
2003), which suggests relatively late reionization; detailed model- 
ing of the proximity zones around z > 6 quasars fWvithe & Loebl 
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2004bl: lMesinger & Haimarl2004IWvithe et alj20o3l . which sug- 
gests xm >0.1 at z ~ 6.3; and the lack of evolution in the 
lumi nosity function of Lyq-se lected galaxies over z ~ 5.7- 
6.5 l Malhotra & Rhoads 2004), which suggests reionization at 
an earlier epoch (although this is subjec t to uncertainties about 
the clustering of hig h-redshift galaxies; iFurlanetto et alj [2004c: 
IWvithe & Loebl20ojh . 

Weaving these individual threads into a consistent picture re- 
quires a rather complex history. In the simplest models, reioniza- 
tion proceeds extremely rapidly be cause galaxy formation is nearly 
exponential at high redshifts (e.g., Barkana & Loeb 2001 and ref- 
erences therein). This is incompatible with the combination of 
SDSS and WMAP measurements (at least if we ignore the large 
errors on the latter measurement). The most popular solution is 
to invoke some sort of feedback mechanism to de crease the ioniz- 
ing efficiency of galaxies during reionization (e.g. JWvifhe & Loebl 
2003albl lCeJ2003JbtlHaiman & Holdej2003l) . Feedback may di- 
rectly affect collapsed objects (for example, through enrichment 
by heavy elements) or the IGM itself (for example, through pho- 
toheating). Unfortunately, the underlying mechanism - or mech- 
anisms - remain mysterious, as do their effects on the reion- 
izatio n history (t hough they are unlikely to introduce sharp fea- 
tures; [Furlanetto & Loeb 2005). Obviously, the physics of reion- 
ization, the first luminous sources, and the IGM are deeply in- 
tertwined. Disentangling them requires new methods to study the 
high-redshift universe. 

Fortunately, the existing quasar data hints at one promis- 
ing direction: the two best-studied lines of sight c ontain remark- 
able differences iBecker et all 1200 ll : IWhite et alj 120031) . SDSS 
J1030+0524, at z — 6.28, shows completely saturated absorption 
from z w 5.97-6.18 in all the Lyman transitions. This places a 
lower limit on the (mean) Lya optical depth r a > 9.9 (at the 2a 
level). In contrast, SDSS Jl 148+5251, at z = 6.42, shows a sin- 
gle transmission feature in Lya at z — 6.08, a number of obvious 
transmission gaps in hy/3 over z = 5.95-6.18, and a faint con- 
tinuum in Lya, Ly/3, and Ly7 from z — 6.18-6.33. While at first 
attribu ted t o contamination by a z = 5 interloper. lOh & Furlanetto! 
120051) and lwhite et alj 120051) recently demonstrated that both the 
continuum and t he spikes a re true transmission from the z ~ 6 
IGM. lQh & Furlanettol 1200*51) used the Ly7 trough to place an up- 
per bound on the effective Lya optical depth at z — 6.2-6.3 of 
r a < 15.4 (2a), with the likely value in the range r a ~ 7-11 (the 
large uncertainty comes primarily from the unknown IGM density 
distribution). At lower redshifts, the clear transmission spikes to- 
ward SDSS Jl 148+5251 imply even more significant differences in 
the IGM transmission properties between these lines of sight. Both 
therefore show large jumps in the ionizing background at z ~ 6 but 
are otherwise entirely different. Most interestingly, the differences 
are not "one-shot" affairs that can be localized to a small portion 
of the spectrum: instead they cohere over ~ 100 comoving Mpc 
scales. How can we explain modulation in the transmission over 
such large scales, on which density fluctuations should be tiny at 
z ~ 6? And how can we reconcile it with the sudden and rapid 
evolution of the mean transmission along each line of sight? 

To answer these questions, we must understand how the spa- 
tial distribution of ionized gas evo l ves th rough reionization. Early 
models, such as|Arons & Wingert 1 1972), followed t he growth of 
HII regions around individual galaxies. IFurlanetto et alj {2004a, 
hereafter FZH04) showed that the sizes are actually determined 
by clustering of the underlying galaxy population. Thus studying 
the inhomogeneity inherent to the reionization process can open 
new windows on the physics of reionization. FZH04 showed that 



the ionized bubbles are large throughout reionization and that they 
have a characteristic size. This offers a first step in explaining how 
large-scale fluctuations develop. We review their model in S|2] In 
this paper, we focus on how the HII regions evolve during reioniza- 
tion and especially near its completion. In ij3| we examine how the 
bubbles actually grow and merge with each other in order to eluci- 
date the meanings of jumps in the ionizing background. We show 
that "major mergers" dominate the growth process. As a result, the 
radiation background at any point in the IGM evolves through a 
series of discrete and rapid jumps. However, unlike in previous the- 
oretical treatments (e.g., Gnedin 2000), these jumps do not corre- 
spond to global overlap; rather, they are simply the byproduct of 
inhomogeneous reionization. 

In the original model of FZH04, the bubbles merge until they 
form a single infinitely large object at the completion of reioniza- 
tion. However, by ignoring IGM substructure the FZH04 model 
treats recombinations extremely crudely. iMiralda-Escude et alj 
(2000, hereafter MHR00) presented a complementary model for 
inhomogeneous reionization that focuses precisely on the question 
of recombinations. They argued that locally reionization proceeds 
from low to high-density regions. Once ionized, voids remain so 
because of their long recombination times, but dense clumps (cor- 
responding to filaments or bound structures) rapidly return to neu- 
trality. Thus, within an HII region, the ionized fraction traces the 
fluctuating density field, with the neutral gas confined to smaller 
and denser clumps as reionization proceeds. We describe this model 
in more detail in [g] and we examine some of the subtleties of the 
MHR00 treatment in an Appendix. 

In a sense, FZH04 and MHR00 were designed to address 
different regimes. The former applies best during the pre-overlap 
stage, when the discrete bubble structure from source clustering is 
most important, while the latter applies best during the post-overlap 
stage, when the ionizing background is nearly uniform and the local 
ionized fraction traces IGM density fluctuations. Of course, toward 
the tail end of reionization, both contributions are equally impor- 
tant. We thus show how to integrate the FZH04 and MHR00 mod- 
els in SJ5] The large ionized bubbles of the FZH04 model are of 
course not fully ionized; rather, they contain dense neutral clumps 
that do not affect the bubble topology so long as the mean free path 
between clumps exceeds the bubble radius. As a bubble grows, the 
ionized gas must therefore extend deeper into the dense clumps, 
and the recombination rate increases. Eventually, these internal re- 
combinations consume all the ionizing photons, and the bubble 
stops growing. For the bubble interior, this point is effectively the 
end of "overlap" because bubble mergers have no further effect 
on the ionizing background. Thus overlap is actually a local phe- 
nomenon. In f|6| we show that the resulting saturation has a signif- 
icant impact on the Lyman-series transmission along quasar lines 
of sight toward the tail end of reionization. Finally, we conclude in 

In our numerical calculations, we assume a cosmology with 
n m = 0.3, fi A = 0.7, n b = 0.046, H = lOO/ikms" 1 Mpc" 1 
(with h — 0.7), n = 1, and as = 0.9, consistent with the most 
recent measurements ISpereel et all2003l) . Unless otherwise spec- 
ified, we use comoving units for all distances. 



1 This stands in contrast to the large scale behavior, which traces clusters 
of sources inside overdense regions. 
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Figure 1. (a): Probability distribution of mean free paths in ionized regions. The different curves correspond to different Sj. All curves take f = 12 and range 
from 2 = 8.8-8; we do not take recombinations into account, (b): Evolution of characteristic bubble size (or mean free path) with Xi at a fixed redshift (which 
corresponds to varying Q. The solid, long-dashed, and short-dashed curves are for 2 = 6,9, and 15, respectively. The thick lines assume the FZH04 model; 
the thin lines (shown only for z = 6 and 9) incorporate recombinations as described in i|5] 



2 THE BUBBLE SIZE DISTRIBUTION WITHOUT 
RECOMBINATIONS 



Recent numerical simu lations (e.g., ICiardi et alJ 120031 : 
ISokasian et alj 1200 3l l2004t) show that reionization proceeds 
"inside-out" from high density clusters of sources to voids (at 
least in the standard scenario where galaxies ionize the IGM). We 
therefore associate HII regions with large-scale overdensities. We 
assume that a galaxy of mass m ga i can ionize a mass Cm ga i, where 
( is a constant that depends on (among other things) the efficiency 
of ionizing photon production, the escape fraction of these photons 
from the host galaxy, the star formation efficiency, and the mean 
number of recombinations. Each of these quantities is highly 
uncertain (e.g., massive Popul ation III stars can dramatically 
increase the ionizing efficiency; iBromm et al-lEoOll) . but at least 
to a rough approximation they can be collapsed into this single 
efficiency factor. The criterion for a region to be ionized by the 
galaxies contained inside it is then / co ii > where / co u is 

the fraction of mass bound in halos above some m m i n . We will 
assume that this minimum mass corresponds to a virial temperature 
Tvir = 10 4 K, at which point hydrogen line co oling becomes 
efficient. In the extended Press-Schechter model (Lacey_& Cole 
Il993l) . this places a condition on the mean overdensity within a 
region of mass m, 



2 

mill 



2/ \ll/2 
T (m)] , 



(1) 



5 m > S x (m,z) = S c (z) - V2K(()[a 

where K(Q = erf _1 (l — C _1 )> o- 2 (m) is the variance of density 
fluctuations on the scale m, cr m i n = cr 2 (m m i n ), and S c (z) is the 
critical density for collapse. The global ionized fraction is Xi = 
C/coii,g, where /coll.g is the global collapse fraction. 

FZH04 showed how to construct the mass function of HII 
region s from S x in an analogous way to t he halo mass function 
IPress & Schechterll 19741 bond et alJ ll991). The barrier in equa- 
tion Q is well approximated by a linear function in a 2 , 5 X « 
B(m,z) = Bq + B\o (m). In that case the mass function of 



ionized bubbles (i.e., the comoving number density of HII regions 
with masses in the range mi dm/2) has an analytic solution l Sheth 
fl99aiMcOuinn etall20051: 



m,z=\- 



d In a 



dlnm 



B (z) 



■ exp 



B 2 (m,z) 



2cr 2 (m) 



(2) 



where p is the mean density of the universe. FZH04 showed some 
examples of how the bubble sizes evolve throughout the early and 
middle stages of reionization (see also Fig.|4ji). Figure^ plots the 
distributions near the end of reionization (the era with which we 
will be most concerned) assuming £ = 12. We have phrased the 
size distribution in terms of the mean free path of ionizing pho- 
tons by equating that quantity to the bubble radius: this implic- 
itly neglects absorption within the bubbles, whose importance we 
consider in iJ5] below. Panel a shows the probability distribution 
P(R)dR of the mean free path, or the fraction of ionized matter 
within bubbles of comoving size R ± dR/2. 

Figure illustrates a number of key points about the model. 
First, the bubbles are large: tens of comoving Mpc during the fi- 
nal stages. This is because high-redshift galaxies cluster strongly, 
so small mean overdensities translate to much larger overdensities 
of galaxies. Furthermore, the bubbles attain a well-defined charac- 
teristic size -Rchar at any point during reionization. This occurs be- 
cause of the Gaussian nature of the underlying linear density fluc- 
tuations and because the barrier 5 X is a (decreasing) function of 
m. Near the completion of reionization, over- and underdense re- 
gions sit in large and small bubbles, respectively. In contrast, the 
barrier used in constructing the halo mass function, S c (z), is in- 
dependent of mass, which yields the usual power law behavior at 
small masses. 

The thick curves in Figure^ show how this characteristic 
size depends on Xi and redshift. We vary C, to fix Xi at z = 6, 9, 
and 15 (solid, long-dashed, and short-dashed curves, respectively). 
Regardless of redshift, _R c hai- surpasses 1 Mpc relatively early in 
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reionization (note that this is already much larger than an HII re- 
gion around an individual galaxy) and then increases at an accel- 
erating rate. Bubbles become larger if reionization happens earlier, 
because the underlying galaxy population is more biased at high 
redshifts. However, by the late stages (xi > 0.75), the differences 
are relatively small. For our purposes, the crucial result is that R c h& T 
rapidly approaches infinity when Xi > 0.9. 

As emphasized by FZH04, this simple analytic model makes a 
number of assumptions about the reionization process. But the most 
glaring is that the bubbles are roughly spherical: clearly the inho- 
mogeneous distribution of sources and the cosmic web will intro- 
duce asymmetries throughout reionization. This is not as bad as it 
may sound, because the model actually computes the total mass of 
ionized regions; the radius is therefore only a characteristic value, 
and it may still be relev ant. Simulations like lSokasian et all 1200 3h 
and Ciardi et al. (2003) do seem to produce roughly spherical re- 
gions throughout most of reionization. However, a more serious 
problem is that the topology may reverse itself late in reionization: 
the universe could contains neutral islands embedded in an ionized 
sea, rather than the other way around. This certainly happens in 
simulations, although the timing is suspect because of their small 
sizes. In the FZH04 model, the universe is instead divided into iso- 
lated bubbles separated by thin walls of mostly neutral gas. That is 
clearly an approximation, but it still may be true that i? c h ar gives 
the typical distance to a neutral island. The regime over which the 
model applies must be determined through detailed comparison to 
numerical simulations; however, simulations on the relevant scales 
are not yet available (and the size distribution of HII regions has 
yet to be quantified in any existing simulations!). Another crucial 
assumption is that the recombination rate can be approximated as 
uniform across the entire universe: the purpose of Sj5|is to test the 
validity of this assumption. Along the way we will see that it also 
alleviates the topological problem, and that i? c h ar describes the typ- 
ical mean free path of an ionizing photon. 



3 BUBBLE GROWTH IN THE EXCURSION SET MODEL 

Figure 0; shows that bubbles grow rapidly, especially in the late 
stages of reionization. In this section we will study the mechanisms 
driving that evolution. 



3.1 The progenitor size distribution 

One advantage of the excursion set formulation of the dark matter 
hal o mass function is that conditional mass functions follow eas- 
ily iLacev & Colell993h . The same holds true for the bubble mass 
function. Suppose we consider a region of the universe with over- 
density 8m on scale Mb corresponding to o\. The derivation of 
the mass function follows exactly the same procedure as in FZH04, 
except that the random walk begins at the point (5m, a 'b) so me 
(linear) absorbing barrier has a different amplitude. In particular, 
we can imagine that our initial point corresponds to a bubble of 
size Mb at Zb- in that case, 8m = B(Mb, z&), and the conditional 
mass function represents the bubble progenitor distribution at some 
z > Zb. It is easy to show that 



n b (m, z\M b , z b ) 



P_ 

ty m 2 

2 [B(M, 



d In a 



dlnm 



b,Z) 



B(M b ,Zb)] 



Note that this has the same structure as the unconditional rib(m); 
the only differences are the substitutions 



Bo(m, z) 
B(m, z) 



B(M b ,z)- B(M b ,z b ) 
B(m,z)- B(M b ,z b ) 



(4) 



2 2 



which simply represent translating the origin in the diffusion prob- 
lem. It is also structur ally similar to the halo progenitor distribution 
lLacev & CoIelll993h . with the exception of a linearly-increasing 
barrier, of course. 

Figure [2] shows some example distributions at z — 13 as- 
suming £ = 40. Each panel shows the fraction of the universe in 
bubbles of a given radius, normalized by the mean ionized frac- 
tion Xi = 0.43. The solid curve shows the unconditional n b (m, z). 
The dotted, dot-dashed, long-dashed, and short-dashed curves show 
n b {m,z\M b ,Zb) if z b = 12 and M b = 10 13 , 10 14 , 10 16 , and 
10 18 M ; these correspond to R b = 3.9, 8.4, 38.9 and 180 Mpc. 
Panels a and b assume Zi = 12 (when Xi = 0.68) and Zi = 11.2 
(Xi = 0.96), respectively. Note that the areas under these curves 
are not normalized to unity, because the local ionized fraction at 
z = 13 is greater than the mean value. There are a number of 
interesting features to the progenitor distribution. First, large ion- 
ized regions tend to contain larger-than-average bubbles at earlier 
times. (This is not true for small bubbles, because progenitors can- 
not be larger than the final object.) This is of course because such 
bubbles encompass overdense regions where structure formation 
begins earlier than average, so "local overlap" occurs at an ear- 
lier time. However, one might naively expect this bias to become 
ever stronger as the final bubble size increases. In fact the con- 
ditional rib approaches a constant form, because of — > and 
B(Mb, Zb) —> Bo(zb) on large scales, so equation J5J approaches 
a constant functional form. Physically, mass fluctuations on such 
large scales are extremely small. Thus, all large bubbles have essen- 
tially the same underlying source population and hence essentially 
the same rib - It is only at the end of reionization, when density fluc- 
tuations appear amplified by the proximity to the phase transition 
at full ionization, that such low-amplitude large-scale modes con- 
tribute to differences in the bubble population. By that time even a 
small number of extra sources can connect the already-large bub- 
bles into one fell beast. Earlier on, this tiny bias does not manifest 
itself as clearly because the extra sources only ionize a small phys- 
ical volume. 

As a consequence, the progenitor distribution remains rela- 
tively unbiased if Zb is near full reionization. Mathematically, this 
occurs because B(Mb, Zb) ~ 0. Physically, the mean overdensity 
in ionized regions is extremely small and the boost in the num- 
ber density of ionizing sources is modest - even if the linear bias 
is large; thus their progenitors are nearly average parts of the uni- 
verse. This can also be seen in FigureQ which shows how the size 
distribution becomes narrower as Xi — > 1. Thus, the progenitor bias 
shows up most strongly early in reionization. These conclusions are 
independent of our choice for f . 

Another way of looking at the progenitor bias is to consider 
the evolution of the ionized fraction within bubbles, which is sim- 
ply proportional to their collapse fraction. In the linear barrier ap- 
proximation, 



(<r> 



,2\3/2 



i(z\M b , z b ) = C erfc 



8 c (z)- B{M b ,Zb) 



V 2 Kin-- 6 2 ) 



(5) 
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Figure 2. Progenitor bubble distributions, (a): The solid curve shows the 
overall size distribution at z = 13 if £ = 40 (corresponding to Sj = 
0.43). The other curves show the size distribution at z = 13 inside an 
ionized bubble that at z b = 12 (& = 0.68) has M b = 10 13 , 10 14 , 10 16 , 
and 10 18 (dotted, dot-dashed, long-dashed, and short-dashed curves, 
respectively), (b): Same, except 2% = 11.2 (a;, = 0.96). 



Figure 3. History of the ionized fraction within bubbles relative to the 
global mean. The upper (thick) and lower (thin) sets of curves take z b = 12 
and 11.2; the ordinate R refers to the final bubble size at 2%. Dot-dashed, 
dashed, and solid curves are for z = 18, 15, and 13. The thick and thin 
dotted curves show the normalized bubble size distributions at 2% = 12 and 
11.2, respectively. All assume £ = 40. 



3.2 Merger rates 



The upper sets of curves in Figure[3]show Xi(z\Mb, Zb)/xi(z) at 
several redshifts if £ = 40 and Zi = 12 (thick curves) and 11.2 
(thin curves). The dot-dashed, dashed, and solid curves are for z — 
18, 15, and 13, respectively; these have Xi = 0.034, 0.17, and 
0.43. For reference, the thick and thin dotted curves in Figure [5] 
show the bubble size distributions at zi = 12 and 11.2. Obviously, 
large bubbles tend to be ionized earlier than average, while small 
bubbles tend to be ionized later. The latter occurs because bubbles 
much smaller than the characteristic size appear only in underdense 
regions, which contain few ionizing sources. At large masses, the 
ratio approaches a constant. This is clear from equation because 
B(Mb) — > -Bo and Ob <C m m for Mb — » oo. Physically, the 
explanation is the same as for the progenitor distribution: large- 
scale modes only add a small number of sources that cannot ionize 
much of the IGM but can connect existing HII regions. Thus these 
modes only affect the topology when Xi ~ 1. For the same reason, 
we also see that the bias increases if we select bubbles early on: 
when Xi ~ 1, ionized regions have nearly the mean density and 
hence their histories are near average as well. On the other hand, 
for a given bubble the ionization bias is largest early in reionization 
when galaxies are the most biased, because high-redshift galaxies 
are so far off on the nonlinear-mass tail that a small boost to the 
argument of the error function increases x% by a large amount. 

We have assumed throughout that cr m i n is fixed at a mass scale 
corresponding to T v \ r — 10 4 K and that £ is independent of mass. 
In reality, a variety of feedback effects may decrease the star for- 
mation efficiency in shallow potential wells, effectively increasing 
m-min or even turning £ into an increasing function of mass. This 
would only exacerbate the ionization bias discussed above, so long 
as we normalize to a fixed Xi(z), because it would increase the 
mean bias of the sources. 



We can use the above expressions to derive approximate merger 
rates for ionized bub bles. The procedure is exactly the same as in 
lLacev&Coielll993l) : the result is 



d 2 p(mi, niT, t) 
drri2 dt 



2 1 
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dlnB 
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■/a?)3/2 



1 



1 



(6) 



where d 2 p(mi,mT, £)/dm2 dt is the probability that a given halo 
of mass mi will merge with a halo of mass mi = mr — mi per 
unit time. We will also define a merger kernel 



Q(mi, rri2,t) 



d 2 p(mi,iriT, t) 



n b (m2,t) dm,2 dt 



(7) 



so that the total rate at which bubbles of mass m\ merge with those 
of mass mi per unit volume is n b (mi)nb(rri2)Q{rn\,m2) ■ 

Before proceeding further, we note that this expression suffers 
from the same li mitations as the usua l extended Press-Schechter 
merger rates (see iBenson et alj (2005 for a detailed discussion). 
The most obvious shortcoming is that Q{m\,m2) is not sym- 
metric in its arguments, which becomes especially important for 
large mass ratio mergers. Moreover, the merger rates do not prop- 
erly conserve mass. At root, the difficulty emerges because the 
excursion set formalism considers the trajectories of independent 
points in the IGM rather than grouping them into discrete halos 
which can then merge in a well-defined manner. The mathemati- 
cally self-consistent solution is to invert the Smoluchowski equa- 
tion, but that is numerically challenging. Fortunately, despite their 
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Figure 4. (a): Size distribution of bubbles at z = 12 if C = 10, 20, 30, 40, 
and 50 (dotted, dot-dashed, short-dashed, long-dashed, and solid curves, re- 
spectively), (b): Merger rates at z = 12 for mi = 10 14 Mq. The different 
curves correspond to the same choices for f as in the top panel. The ordinate 
is defined in eq. j8j. 



inconsistencies extended Press-Schechter halo merger rates do of- 
fer useful insights into galaxy formation. We hope that the bubble 
merger rate given by equation JSJ does the same. To correct for 
the asymmetric merger kernel, we will use the symmetrized ver- 
sion 2Q aym (mi, mz) = Q(mi,m 2 ) + Q(m 2 , mi). Note that the 
aforementioned difficulties may be less severe for ionized bubbles 
because of their small-mass cutoff. 

Figure [4] shows some sample merger rates at z = 12. 
The top panel shows the bubble size distribution if £ = 
10, 20, 30, 40, and 50 (from left to right). These correspond to 
Xi = 0.17, 0.34, 0.51, 0.68, and 0.85. The bottom panel shows 
the merger rate for bubbles with mi = 10 14 Mq, The ordinate is 

V(m 2 ) 



vf ^ ldV 
V( mi ) — 



, m 2 rif,(m 2 ,; 
V (mi) 

x Q B ym{rn 1 ,m 2 ,t) 



(8) 



which gives the fractional volume accretion rate by which bubbles 
of mass mi grow through mergers with bubbles of mass m 2 (in 
logarithmic mass intervals). Clearly, regardless of Xi, ionized re- 
gions tend to merge most rapidly with bubbles slightly above the 
current characteristic radius, so mergers with larger systems dom- 
inate the accreted mass regardless of the bubble of interest. The 
total mass accretion rate increases rapidly throughout reionization 
because there is more gas in ionized bubbles at late times and be- 
cause each bubble is larger. Thus HII regions grow at an acceler- 
ating rate as Xi — > 1 (see also Fig.Q>). Finally, note that the frac- 
tional accretion rate approaches a fixed asymptotic form for small 
bubbles. This is because equation |6j becomes independent of m 2 
(jriT ~ mi and o\ ~ a 2 ) when m 2 <C mi. 

Figure [5] shows similar results, except we fix £ and vary 
mi. We show the volume accretion rate (in dimensional co- 
moving units), again as defined in equation {5}. The solid, 
dot-dashed, short-dashed, and long-dashed curves take mi = 



Figure 5. Merger rates at z = 12 for mi = 10 11 , 10 12 , 10 14 , and 
10 16 Mq (solid, dot-dashed, short-dashed, and long-dashed curves, respec- 
tively). The ordinate is defined in eq. JS). {a): C = 20. (b): C = 40. 



10 , 10 , 10 14 , and 10 16 Mq, respectively. These have Ri = 
0.84, 1.8, 8.4, and 38.9 comoving Mpc. The two panels assume 
C = 20 (Xi = 0.34) and £ = 40 (i, = 0.68); Figure^ shows the 
corresponding n(,(m). Interestingly, large bubbles grow faster, in 
absolute terms, than small bubbles. The disparity is greatest for the 
accretion of small bubbles by large ionized regions: apparently, as 
these expand they swallow up small bubbles quite rapidly. On the 
other hand, the rate at which a given bubble merges with large HII 
regions has a much weaker mass dependence. In particular, if we 
had normalized these curves by Ffmi)" 1 (as in Fig.|4j, we would 
find that small bubbles grow fractionally at much faster rates than 
large bubbles, because both tend to merge with ionized regions of 

Size -Rchar- 

These merger rates have interesting implications for our un- 
derstanding of the overlap process. Because i? c har is well-defined, 
mergers with objects of this size dominate bubble growth - re- 
gardless of the bubble under consideration. This implies that small 
bubbles, which grow in voids, do not merge with each other: in- 
stead they overlap only when a large bubble, from a distant clus- 
ter of sources, sweeps past them. Such behavior is consistent with 
the "inside-out" nature of reionization. It also helps to explain the 
rapidly increasing R! c har seen in Figure^: most of the universe lies 
inside bubbles near -R c har, which merge with bubbles of slightly 
greater sizes at ever-increasing rates. 

Most importantly, this picture also affects our interpretation 
of the evolving radiation background. Figure|6|illustrates how. We 
have generated an ensemble of random walks in the (a 2 , 8) plane 
and tracked their bubble histories. We fixed Xi = 1 at z = 6, as- 
sumed constant £, and compared each random walk to our FZH04 
barrier at intervals of Az — 0.05. We begin at z = 13, when 
Xi — 0.065. The solid and dashed curves are for models with 
and without recombinations, respectively (to be discussed in J5]be- 
low); for now the middle stages are most relevant. The dotted line 
shows -R c har for reference. Even though Xi grows smoothly in this 
model, the process is extremely discontinuous if viewed from any 
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Figure 6. Bubble histories for several randomly generated trajectories. The 
abscissa shows the bubble radius surrounding each point as a function of 
Xi and z; we fix £ so that reionization completes at z = 6. The solid lines 
include the recombination limit (see J5]below) while the dashed curves use 
the unmodified FZH04 model (the two are coincident until Sj > 0.9). The 
dotted and dot-dashed lines show -R c har and iimax, respectively. 



given point: although the local emissivity grows only slowly, the 
point joins progressively larger bubbles as time goes on. Moreover, 
large bubbles always dominate the volume accretion rate, so each 
merger is "major" (i.e., it typically increases the volume by a factor 
of two or much more) and causes a substantial jump in the bub- 
ble radius. Note that this is even true when bubbles contain many 
ionizing sources - the variance does not decrease like a Poisson 
process. Because (neglecting recombinations) the bubble radius is 
also the horizon out to which galaxies are visible, this in turn im- 
plies that the local ionizing background evolves through a series 
of sharp, large amplitude jumps. Obviously we must be careful not 
to interpret these local jumps in the ionizing background as sig- 
natures of global overlap. This contrasts with the stand ard picture 
from simulations (Gnedin 2000; Razoumo v et alJl2002l) . in which 
a jump in the ionizing background was attributed to the comple- 
tion of reionization. These simulations were limited to small boxes 
with sizes ~ 1-10 Mpc. As such, they are much smaller than the 
typical HII regions at the end of reionization, and the simulation 
jumps occurred when a bubble encompassed the entire box. But 
because the small boxes are not characteristic of the entire universe 
(Barkana & Loeb 2004; FZH04), this event need not correspond to 
global overlap. Our models show that interpreting such jumps is 
more subtle. 



4 THE IGM DENSITY DISTRIBUTION 
4.1 The MHR00 model 

To this point, we have essentially ignored recombinations in- 
cludes only a factor accounting for the mean integrated number per 
baryon). In the following sections, we will consider how recombi- 
nations in an inhomogeneous IGM affect the bubble topology. To 



begin, we will review a simple analytic model for the IGM density 
distribution and its ionization state originally due to MHR00. They 
used numerical simulations to measure the volume-weighted den- 
sity distribution of IGM gas, Py(A), where A = p/p, and found 
a good fit at z ~ 2-4 with 

£V(A)dA = AoA^ 



cxp 



(A 



-2/3 



Co) 2 



2(2V3) 2 



dA. 



(9) 



Intuitively, the underlying Gaussian density fluctuations are mod- 
ified through nonlinear void growth and a power law tail at large 
A. MHR00 argued that the form could be extrapolated to higher 
redshifts in the following way. First, 5o essentially represents the 
variance of density fluctuations smoothed on the leans scale for an 
ionized medium; thus <5o oc (1 + z) -1 at high redshifts, as their 
simulation found for z = 2-4. The power-law exponent /3 deter- 
mines the behavior at large densities; for isothermal spheres, it is 
P = 2.5. The remaining constants (Ao and Co) can be set by de- 
manding proper mass and volume normalization. 

The next step is to determine how this gas is ionized. MHR00 
assumed that there exists a density threshold Ai such that all the 
gas with A < Ai is ionized while all gas with A > A; remains 
neutral. This provides a reasonable description of shielding in dense 
regions, if those regions can be considered to be isolated clumps in 
which the density increases inwards. In that case, the radiation field 
will ionize the outskirts of the cloud until r ~ 1. Because of the 
density gradient, this skin corresponds to our threshold Ai. 

Given this simple prescription, the recombination rate of ion- 
ized gas, denoted A (note the change from MHR00 in order to avoid 
confusion with the bubble radius), is 



A(Ai) = A t 



dAPy(A) A = A U C, 



(10) 



where A u is the recombination rate per hydrogen atom in gas at the 
mean density, 



A u = a A (T)n e 



(11) 



and C is the effective clumping factor. Here, following MHR00, we 
use the case-A recombination rate (see also lMiralda-Escudel2003l) . 
a A = 4 x 10~ 13 cm 3 s" 1 (at T = 10 4 K). Case-A is appropriate 
because most recombinations occur in gas with A ~ Ai, which 
presumably lies on the edge of dense neutral clouds. Thus the sec- 
ondary ionizing photons that follow recombinations to the ground 
state will find themselves in the dense, neutral clouds and do not 
ionize the diffuse IGM in which we are interested. In this picture 
C increases rapidly as the universe becomes more and more ion- 
ized; thus, the first phase of reionization is more efficient than later 
stages, when we must clear out dense, rapidly-recombining blobs 
of gas. 

MHR00 also offer a prescription for determining Ai, the mean 
free path of ionizing photons. In their model, it equals the mean 
distance between clumps with A > Ai along a random line of 
sight, which is approximately 



A, 



A [l-iV(Ai 



-2/3 



(12) 



Here Fv(Ai) is the fraction of volume with A < Ai and Ao is 
a (redshift-dependent) normalization factor. Formally, this expres- 
sion is valid only if the number density and shape (though not total 
cross section) of absorbers is independent of Ai. This is obviously 
not true in detail for the cosmic web. However, MHR00 found that 
it provided a good fit to numerical simulations at z = 2-4 if we set 
XoH(z) = 60 km s _1 (in physical units). We will extrapolate the 
same prescription to higher redshifts. 
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While it does a remarkable job matching many aspects of the 
Lya forest at moderate redshifts, applying the MHROO model to 
high redshifts involves a number of potentially dangerous simplifi- 
cations and extrapolations. In the Appendix, we examine the model 
in some more detail and illuminate many of its features. For ex- 
ample, we show that equation i\2\ probably overestimates A; by 
up to a factor ~ 2 because of accumulated photoelectric absorp- 
tion by low-column density systems. We urge the reader to keep in 
mind the approximate nature of our model. We emphasize that bet- 
ter measurements with high-resolution simulations will be useful 
in the future. 

4.2 A minihalo model 

One possibly unwarranted assumption of MHROO is the explicit 
smoothing of the IGM on the Jeans scale for T = 10 4 K gas. While 
appropriate well after reionization, this may not apply during that 
epoch. If the IGM remains cool before reionization, small collapsed 
objects called "minih alos" could form; they will eva porate dur- 
ing reionization iBarkana &Loeb 1999; Shapiro et al. 2004) but 
may constitute a substantial photon sink before they are completely 
destroyed (Haiman et al. 2001; Barkana & Loeb 2002; Iliev et al. 
2005). However, the importance of minihalos remains controver- 
sial: their formation relies on fluctuations on extremely small 
scales, and even a modest amo unt of hea ting will strongly sup- 
press their abundance <Oh&Haimanl2 003). For transparency, and 
because minihalo formation is so unconstrained, we therefore use 
an extremely simple model for a universe in which minihalos con- 
trol the recombination rate. We suppose that a fraction / m h of the 
IGM has collapsed into dense blobs with mass M m h and mean den- 
sity A m h ~ 18-7T 2 . We will assume that these minihalos remain 
neutral, randomly distributed sinks for ionizing photons. Then the 
number density of minihalos is n m h = fmhP/M m h, each with ra- 
dius Rmh oc (M m h/A m h) 1//3 . The mean free path of an ionizing 
photon is 

; _ 1 M mh V /3 ^0.05^ 

X (l8^J (,WJ MPC U3) 

in comoving units. The fraction of ionizing photons surviving to a 
distance R from a source is f s = 1 — exp(— 7?/Z m h). 



5 BUBBLES WITH RECOMBINATIONS 
5.1 The recombination barrier 

We are now in a position to incorporate recombinations into the 
FZH04 model. Our fundamental underlying assumption is that the 
bubble distribution depends on the underlying density field: HII 
regions correspond to (large-scale) overdense regions of the uni- 
verse. We will therefore need to write the recombination rate as a 
function of the (smoothed) mean density. Conceptually, we must 
modify equation 1101 in two ways. First, the recombination rate per 
baryon increases with the mean density of the region: A u oc (1+5). 
Second, denser regions can be thought of as sub-universes with a 
larger matter density. In such regions, structure formation proceeds 
more rapidly than average [H(z) oc \ZQ m at high redshifts], so 
they are slightly clumpier than average. This adds a second layer to 
the dumpiness calculation: in principle fV(A) — ► Pv(A\5,R). 



For simplicity, and because the overdensities of bubbles are modest 
during most of reionization (one implication of their large sizes), 
we will neglect the latter effect and assume Pv (A) is independent 
of environment. 

Now suppose that we have a bubble of size R. In order for the 
bubble to grow, ionizing photons must be able to reach its edge. 
This demands that A,: > R, setting the threshold A, to which we 
must have ionized the gas inside the bubble through equation I12i . 
Obviously, larger bubbles contain more highly ionized material, be- 
cause their photons must propagate farther. At the same time, the 
ionized outskirts of dense clumps are recombining, and we must 
counteract these recombinations in order for any photons to reach 
the edge of the HII region. Thus for bubbles to continue growing, 
the rate at which they ionize their constituent gas must exceed the 
rate at which it recombines: 

>A U C'(5,R). (14) 

The condition Aj = R sets A*; with our assumption that Py(A) is 
independent of the mean bubble density, we can then write 

C'(8,R) = C(R)(l + 8). (15) 

Note that bubbles are not completely ionized in this picture, be- 
cause each contains some fraction of dense neutral gas. The quan- 
tity Xi should be thought of as the fraction of matter contained in 
HII regions rather than the fraction of hydrogen atoms that are ac- 
tually ionized. 

Equation 1 141 differs from the FZH04 condition / co ii > 
which requires only that the sources inside produce enough photons 
to ionize every hydrogen atom when integrated over the history of 
the universe (or at best to counteract the cumulative number of re- 
combinations). The condition of equation 1 141 instead requires that 
the sources produce enough ionizing photons to counteract recom- 
binations at that instant. Clearly these are two limiting cases: in re- 
ality the sources must produce enough photons to ionize each atom 
plus cancel ongoing recombinations. It is thus nothing more than 
the canonical Stromgren radius, in which recombinations exactly 
balance ionizations. Separating the problem into these two limits is 
useful so long as one or the other clearly dominates over most of 
parameter space, which as we shall see below is indeed the case. We 
note that MHROO made the same approximation in their Figure 7; 
our treatment simply extends theirs to include the inhomogeneous 
distribution of (clustered) sources throughout the universe. 

Equation J14I can be solved numerically to yield the minimum 
S r as a function of scale M required for ionizing sources to counter- 
act recombinations inside a bubble. In the excursion set formalism, 
a trajectory must then cross St = ma,x(S x , S r ) in order to be part 
of an ionized bubble. We compare these barriers in Figure EJi. The 
three monotonically decreasing curves show 8 X , the FZH04 barrier, 
at Xi = 0.95, 0.82, and 0.49 (solid, dotted, and dashed lines, re- 
spectively). 2 The thick nearly vertical curves show <5 r for the same 
neutral fractions. On small scales, recombinations provide a much 
less stringent limit than simply ionizing the gas: small bubbles have 
a low density threshold A,, because the photons need not propagate 
far. Thus the recombination time inside each bubble is long. How- 
ever, as we increase the required A;, A increases more rapidly than 
the collapse fraction (see, e.g., Fig. 2 of MHROO) and 5 r becomes 
the limiting factor. The steepness of 8 r is a result of the power-law 
form of equation {5J at large A, which assumes that dense clumps 

2 As shown in FZH04, these barriers are nearly linear when expressed as 
functions of c 2 (m). 
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Figure 7. Barrier height as a function of comoving radius. The solid, dotted, 
and dashed curves have Xj = 0.95, 0.82, and 0.49, respectively, at z = 6. 
Within each set, the thin curves show S x (from FZH04) and the thick curves 
show <5 r (from this paper), (a): Uses MHR00 to compute S r . (b): Uses the 
minihalo model described in the text to compute 8 r . 

are isothermal spheres, p oc r~ 2 . Increasing A; requires a substan- 
tial increase in Aj and hence in the recombination rate. Meanwhile, 
increasing d/ co n / dz is difficult, so the required underlying density 
grows rapidly. The crossover between the two limits depends on 
the ionized fraction, but it typically occurs at R > 10 Mpc. FigureQ 
shows that the bubbles do not approach this scale until Xi > 0.75, 
so the recombination limit only affects the bubble size distribution 
late in reionization. : explanation of the title We thus see that re- 
combinations play an increasingly crucial role as the bubbles grow 
larger and larger. They constitute a progressive tax on the ionizing 
sources, with richer clusters paying more because their bubbles are 
intrinsically larger. 

Figure Q suggests that the crossover scale i? max , where 8 X — 
St, does evolve as reionization proceeds. However, the rate of evo- 
lution is quite modest. We show this explicitly in Figure^. The 
thin curves (shown for z — 6 and 9) plot the maximum of 7? c har 
and J?max; the turnover at Xi>0.9 is due to the recombination 
limit. Thus, HII regions continue to grow after recombinations be- 
come important, but only slowly - the evolution toward Xi ~ 1 is 
completely different from a case with no recombinations, in which 
the bubble growth is unbounded. This slow growth is similar to 
what we expect in the post-overlap universe, because bubble merg- 
ers no longer increase the ionizing background (see H5.3t , 

We can apply similar arguments to a universe in which miniha- 
los control the mean free path. In this case, the fraction of photons 
propagating to the edge of the HII region is f s (recall we assume 
that any photon impacting a minihalo is lost irrevocably). Thus, the 
efffective dumpiness is 

C mh (R) = (1 - ,f mh ) 2 exp(R/l mh ), (16) 

where the first factor describes the fraction of mass remaining in 
the diffuse IGM. We again take C' mh (5,R) = (1 + 8)C mh (R), 
which assumes that minihalo formation is unbiased. That is cer- 
tainly not true, but given the simplicity of our treatment it suffices. 



Figure 8. Barrier height as a function of comoving radius. The solid, 
dashed, and dotted curves in each panel are for z = 6, 9, and 12; the mono- 
tonically decreasing curves show 8 X , while the others are the recombination 
barriers S r (assuming the MHR00 model), (a): Xi = 0.9. The dot-dashed 
curves set XoH(z) = 30 km s _1 for S r (see text), (b): Xi = 0.5. 

(Note also that I liev et al 1200 5 argue that the effective bias of mini- 
halos around a source is more or less independent of scale, which 
implies that we are not ignoring any crucial scale dependence.) Fig- 
ure 0? shows 5 r for our minihalo model with / m h = 0.05 and 
A^mh = 10 6 Mq, We see that these parameters provide a more 
stringent constraint on the bubble sizes than MHR00: the maximum 
size is limited to ~ 1 Mpc during the middle stages of reionization 
and somewhat under ~ 10 Mpc toward the end. Thus, in principle, 
minihalos can strongly affect the topology of reionization; however, 
at least in this simple treatment, they do not completely destroy the 
large-scale fluctuations. Clearly, measurements of the distribution 
of ionized gas can provide powerful limits on their formation. 

Figure|S|shows the barriers for z = 6, 9 and 12 (solid, dashed, 
and dotted curves, respectively) at fixed Xi = 0.5 and 0.9. As 
pointed out by FZH04, the simple S x barrier changes relatively lit- 
tle with redshift so long as we fix Xi\ the only difference is the slow 
decline of m m i n with cosmic time. However, recombinations do 
evolve strongly with redshift, because they depend directly on the 
dumpiness. Interestingly, the MHR00 5 r not only crosses S x at dif- 
ferent scales as a function of z, but it also changes slope. Although 
the parametric form of Pv (A) remains constant with redshift, S r 
is sensitive to a different portion of JV(A) at each redshift. Be- 
cause the dumpiness decreases at high redshifts, the decline in the 
net recombination rate actually allows the bubbles to grow larger, 
even though A u itself increases. In any case, we note that if the 
MHR00 model is accurate, recombinations will only limit the bub- 
bles on tens of Mpc (or even 100 Mpc) scales: with this model, 
reionization is inevitably accompanied by large-scale fluctuations 
in the ionizing background. 

As described in the Appendix, the MHR00 prescription may 
overestimate the mean free path of ionizing photons by up to a 
factor of two. The dot-dashed curves in Figure |8jj show that set- 
ting \qH(z) — 30kms _1 would, not surprisingly, decrease the 
maximum bubble size by about the same factor (they are otherwise 
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constructed in the same way as the usual S r curves). Because of 
the uncertainties inherent to existing models of the IGM density 
field we cannot quantitatively determine the maximum bubble size 
with much confidence; however, eliminating fluctuations on scales 
> 10 Mpc requires doing a fair amount of violence to the model, 
such as adding substantial numbers of minihalos. 

The scale at which S r = 5 X is also quite sensitive to the other 
input parameters. This is a consequence of its steepness: order unity 
changes to the right-hand side of equation J14I require correspond- 
ing changes to d/ co ii/d;z, which (due to its weak scale dependence) 
cannot be done without significantly altering the bubble size. For 
example, in equation Jl H . we have assumed that case-A recombi- 
nation is appropriate, because the re-emitted photons will generally 
be lost inside dense neutral regions. If we instead assume case-B 
recombination, A U (B) ~ 0.6^(^4). We find that this imposes a 
limiting size about twice as large as the case-A values at z = 6. 
This is another reason why numerical simulations will eventually 
be required to fully understand reionization. 



5.2 The bubble size distribution with recombinations 

We have argued that treating the excursion set barrier as the max- 
imum of S x and 8 r is a good approximation for the effects of re- 
combinations on bubble growth. Armed with this barrier, we can 
now consider its effects on the bubble size distribution. From Fig- 
ures0and|8| it is obvious that recombinations are unimportant until 
some scale i? max is reached, where S x = S r and above which S r 
rises rapidly. The slope in this regime depends on redshift, but it 
is always much steeper than the nearly flat 8 X barrier. Thus, the 
primary effect of recombinations will be to prevent bubbles grow- 
ing larger than i? max . Early in reionization, when i? max 2> -Rchar, 
recombinations can be safely ignored. Later on, we expect the bub- 
bles to saturate near i? max . The trajectories in Figure |6| bear out 
this intuition. The dot-dashed curve shows i? max (,z). The solid and 
dashed lines show the bubble histories with and without recombi- 
nations, respectively. Early on, the two sets of curves are coinci- 
dent because recombinations have no effect on the bubbles. How- 
ever, when ii>0.9, R c h&r > fi max and the trajectories depart 
from one another. Without recombinations, the bubbles continue to 
grow; with them, they saturate and cluster around i? max . 

As in FZH04, we must approximate the barrier with a sim- 
ple form in order to estimate n lec (m), the bubble size distribution 
with recombinations. The most straightforward choice is to assume 
that S r is simply a vertical line that crosses S x at 7? max . This is 
actually an excellent approximation in the (a 2 , 8) plane, because 
a is already so small at the relevant scales. We can then solve for 
n lec (m) by noting that there is no barrier at all for R > -R max . The 
probability distribution at this scale will then be a simple gaussian: 



\Rn 



i 



exp 



da, 



(17) 



where <r max = a(R laBix ). Our barrier is a vertical line beginning 
at B(i? max ), so any trajectory that lies above the FZH04 barrier 
(when evaluated on the saturation scale) will be incorporated into 
a bubble with R — J? max . The number density of such bubbles is 
then 



^Vrec (?^max) — 



2m n 



erfc 



B(Rn 



(18) 



(note that this is a true number density with units inverse comoving 
volume). 




10 10 2 
R (comoving Mpc) 

Figure 9. Bubble size distributions at z = 8 for the unmodified FZH04 
model (thin curves) and including recombinations (thick peaked curves). 
The dotted, short-dashed, long-dashed, and solid curves have Xi = 
0.51, 0.68, 0.84, and 0.92, respectively. The filled hexagons show the 
fraction of the ionized volume in bubbles with R = i? max for the model 
with recombinations. 



Trajectories with <5(i? max ) < B(ii max ) continue their ran- 
dom walks until they cross the FZH04 barrier on smaller scales. 
The result must be independent of the trajectory at R > 7? max ; 
we only care about its value where the barrier begins. The resulting 
mass function is [excluding the true recombination-limited bubbles 
ofeq. {18}]: 

>-B(i? max ) 

rw(m,z)= / d(5p(<5j.R max ) n b (m,z\8, m max ,z),(19) 



where ri),(m, z\S, m max , z) is the conditional mass function for a 
trajectory that begins its random walk at the point (o"n lax , 5). This 
function is identical to equation J5} except with the replacement 
B(A'It, Zb) — * S because the trajectory starts at an arbitrary den- 
sity rather than in a specified bubble. In other words, the net mass 
function is the weighted average of the conditional mass functions 
evaluated over all densities smoothed on the scale 7? max . 

Although the integral in equation I19i can be solved ana- 
lytically, the result is complicated and far from illuminating. We 
therefore simply show some example size distributions in Fig- 
ure [5] The thick and thin curves show the sizes with and with- 
out recombinations, respectively. We take 2 = 8 and show results 
for x t = 0.51, 0.68, 0.84, and 0.92 (dotted, short-dashed, long- 
dashed, and solid curves, respectively). If the characteristic size is 
much smaller than the recombination limit, n ICC « n& except for 
the large bubble tail (e.g., the dotted curve). As _R c h ar — * -R max , 
the bubbles pile up near that scale. Without recombinations, tra- 
jectories in slightly overdense regions tend to cross 5 X on scales 
slightly larger than i? c h ar . With recombinations added, that is no 
longer possible, delaying their first crossing until just after i? max . 
When 7?ch ar 3> -R max , the pileup becomes more significant and 
a large fraction of the universe lies in bubbles at or just below the 
recombination limit. 

The filled hexagons in Figure [5] show the fraction of the ion- 



© 0000 RAS, MNRAS 000, 000-000 



Recombinations and Reionization 1 1 



ized volume in bubbles with 7? = fl max . This is not equal to 
the fraction of the volume in bubbles with R > i? m ax for the 
pure FZH04 model, because some trajectories cross 8 X at large 
scales but then fall below it (i.e., small voids in wide overdensi- 
ties). If the FZH04 barrier were constant, exactly half as many tra- 
jectories would lie in saturated bubbles as lie in large bubbles for 
the pure FZH04 model. 3 With our linear barrier, the difference is 
slightl y larger [comp are our equation 1 1 81 with equation (11) of 
McO uinn et all20 051. The other half of the trajectories that would 
be in large bubbles continue to smaller scales; most fill up the peak 

just below -Rmax. 

For the same reason, the total ionized fraction with recombi- 
nations is actually slightly smaller than C/coii, g : a small fraction of 
trajectories, corresponding to deep voids embedded in large scale 
overdensities, would have crossed the FZH04 barrier only on scales 
R > -Rmax - Fortunately, this is a small correction - most likely less 
important than the residual neutral gas within the bubbles - so we 
will not examine it further here. 



5.3 The "time of overlap" 

Thus, recombinations impose a saturation radius R m ax beyond 
which bubbles cannot grow, as ionizing photons are consumed in- 
ternally rather than reaching the edge of the bubble. From the point 
of view of a single point in the IGM, reionization is complete once 
that point joins such a Stromgren sphere. Even if another HII re- 
gion forms contiguous to its host bubble, the ionizing background 
will change only slightly. "Overlap" effectively occurs at different 
times throughout the IGM and should be considered a local pro- 
cess. After saturation, the mean free path (and hence the ionizing 
background) grows much more slowly (see Figs.Q? and|5J, even 
though the u niverse as a wh ole is not fully ionized. According to 
the picture of lGnedir] feOOOl) . such slow growth is characteristic of 
the post-overlap phase. Instead, we argue that the it is characteristic 
of saturation, which occurs throughout the end of reionization. Our 
model provides a well-defined way to compute this distribution of 
overlap times. We let a; rcc be the fraction of points contained in 
bubbles with R > R m ax (computed using the standard FZH04 bar- 
rier). As shown in the previous section, this is essentially the same 
as the fraction of trajectories in recombination-limited (or nearly 
recombination-limited) bubbles. 

Figure 1101 shows a; rcc as a function of x% at three different 
redshifts [z = 6, 9, and 12 for the thin solid, dashed, and dotted 
curves, respectively) assuming the MHR00 model for the mean free 
path. In all cases, only a small fraction of the IGM reaches satura- 
tion early in reionization: none of these large bubbles form before 
Xi — 0.5, and at least 50% of the universe reaches local overlap in 
the range 0.9 < Xi < 1. This fraction increases at higher redshifts: 
as seen in Figure [8] the lack of substructure in the high-redshift 
universe allows the bubbles to grow much larger at z = 12 than 
at z — 6, so the saturation limit occurs at much larger i? m ax- This 
range in x% corresponds to a small (but far from negligible) redshift 
interval. To fix numbers, the top axis of Figure TlOl shows the redshift 
interval from complete reionization for the z — 6 case, assuming a 
constant efficiency f. In this CilSC. Xrcc ~ 0.5 when Az w 0.3, with 

3 This is the s ame factor of two that bedeviled the original 
iPress & SchechteJ <1974 derivation of the halo mass function. It can be 
most easily understood through the excursion set formalism, because the 
barrier absorbs both those trajectories that continue ab ove it and those that 
reflect off of it IChandrasekhaill943t:lBond et alll99lll . 
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Figure 10. Fraction of the universe in Stromgren spheres; i.e., contained 
in FZH04 bubbles larger than -R ma x ■ Solid, dashed, and dotted curves are 
for z = 6, 9, and 12, respectively. The thin lines use the MHR00 density 
model. The thick solid curve assumes that A; is half the MHR00 value at 
2 = 6. The top axis shows the redshift offset from the time of full ionization 
(if that occurs at z = 6), assuming that f is fixed. 

a long tail to even higher redshifts. The redshift interval is some- 
what smaller at higher redshifts, both because R! max is larger and 
because / co n increases more rapidly (in relative terms). : moved 
up from next paragraph Interestingly, the interval increases if we 
decrease the Stromgren radius: the thick solid curve is the distribu- 
tion in overlap times for z — 6 in a model where we take the mean 
free path of ionizing photons to be one half of the MHR00 value of 
equation 1 12i : according to the Appendix this is a reasonable guess 
for the cumulative effect of low-column density systems. By de- 
creasing Rmax, bubble saturation begins much earlier: a; rcc « 0.5 
at Az « 0.6, with a larger tail. Essentially, the bubbles saturates 
earlier in overdense regions, and underdense regions (which would 
otherwise be engulfed by bubbles propagating outward from fila- 
ments) take longer to catch up. 

Interestingly, our model shows that we cannot identify the 
thickness of the reionization surface with the scale R m ax- Naively, 
one might expect the bubble size to provide a characteristic "radius 
of curvature" for the reionization epoch, with larger bubbles result- 
ing in more cosmic variance in the time of overlap. Iwvithe & Loebl 
(2004a) argued that the width of the overlap surface comes from 
balancing causality with density fluctuations (quantified in a way 
similar to FZH04); essentially they equated R ma x with the width 
of the overlap surface. Our model shows that this is an overes- 
timate because of the spread in overlap times; in fact, a smaller 
Rmax produces greater cosmic variance in the overlap era. Thus, 
the difference in Lyman series transmission between the lines of 
sight to the two highest redshift quasars on ~ 100 Mpc lengthscales 
(Oh & Furlanetto 2005; White et al. 2005) does not necessarily re- 
quire bubbles of size ~ 100 Mpc; instead, it implies that bubble 
growth saturated at considerably smaller length scales; according 
to our model, such variance could come from bubbles < 20 Mpc. 
The strength of our model is that we can compute this variance 
quantitatively. 



© 0000 RAS, MNRAS 000, 000-000 



12 Furlanetto & Oh 




o.oi 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 

0.2 0.4 0_.6 0.8 1 

Figure 11. Recombination rates for inhomogeneous reionization. Each 
curve shows the net recombination rate A scaled to that for a homogeneous 
universe (A u ); the upper thick set includes the bubbles, while the lower 
thin set does not. Solid, dashed, and dotted curves are for z = 6, 9, and 12, 
respectively. The dot-dashed curve includes the variable mean free path at 
z = 9 but not the bubble overdensity. Note that, in the lower redshift cases, 
curves begin at = 2/ co n . 



5.4 The global recombination rate 

MHR00 pointed out that the global recombination rate will in- 
crease throughout reionization as the ionization fronts penetrate 
into denser and denser structures. With their assumption of a uni- 
form ionizing background, they fixed A; by requiring Fv(Ai) = 
Xi (note that we treat Xi as a volume-weighted quantity here). 
Then the global recombination rate A follows directly from equa- 
tion <10t . We show the resulting recombination rates, scaled to that 
of a homogeneous IGM (A u ), by the lower thin set of curves in 
Fieure fTTI The solid, dashed, and dotted curves are for z = 6, 9, 
and 12. In the early or middle stages of reionization, A < A u be- 
cause only the voids are ionized. In this regime, A increases toward 
higher redshifts, because voids are denser (relative to the mean) at 
those times. In the late stages, on the other hand, the dense clumps 
become ionized, A > A u , and A increases toward lower redshifts. 
Because we follow the volume-averaged Xi, we do not enter this 
regime until i; ps 1: dense clumps fill only a small fraction of the 
universe even though they contain a substantial fraction of the mass 
(see Fig. 2 of MHR00). Thus the question of whether reionization 
is most efficient at high or low redshifts depends on the Xi of in- 
terest. At low redshift, the early stages proceed faster, because A u 
is smaller and voids are deeper, but the late stages are much slower 
because dense structures are significantly more common; overall, 
reaching Xi m 1 is easier at early times. 

The shortcoming of this prescription is that it does not include 
the large-scale bias of ionized regions or the full inhomogeneity 
of the reionization process, in which some regions become highly 
ionized long before most of the universe. In our model, the global 
recombination rate is simply a weighted average of the recombina- 
tion rate inside each bubble: 



^4bub = / dm b A(m b ) — n b (m b ) [1 + S x (z)]. (20) 

J P 

Here the last factor is the mean physical overdensity of each bubble. 
We set the recombination rate of each bubble, A(m b ), by equating 
Xi and R b in order to fix Aj. 

In Figure fTTI the upper thick set of curves shows Abub for 
the FZH04 model of equation the redshifts are the same as in 
the lower set. 4 This model increases A over that for uniformly dis- 
tributed sources by a factor of about four, even early in reionization. 
One reason is the biased location of bubbles, which sit in overdense 
regions. Another is that the mean free path (and hence A^) grows 
much faster in big bubbles early on. Because the sources are not dis- 
tributed uniformly, they leave some low-density voids neutral while 
ionizing some higher-density gas, inevitably increasing Ahub- The 
dot-dashed curve illustrates the relative importance of these effects: 
it includes the variable mean free path but ignores the bubble over- 
densities. In that case, Abub is significantly closer to the MHR00 
value early in reionization, when S x (z) is largest, but approaches 
equation 1201 for Xi — + 1 because 8 X — *• in that limit. According 
to Figure ITT1 the two effects conspire to boost the recombination 
rate by a factor of a few throughout reionization. This again il- 
lustrates how recombinations "tax" sources in high-density regions 
more heavily than those in low-density environments by forcing 
clusters of sources to blow large bubbles - with correspondingly 
large recombination rates - even early in reionization. 

We must also note two caveats about this treatment. First, we 
have ignored the recombination limit on the bubble size, and hence 
we overestimate Abub when i? c h ar > i? ma x. In actuality the en- 
hancement will be smaller when Xi > 0.9, but the degree will de- 
pend on the precise value of i? max - Second, if equation I12i over- 
estimates the mean free path, A will increase by a larger factor 
because A, rises. For example, decreasing A; by a factor of two 
increases A by ~ 10-50%, depending on the ionized fraction and 
redshift. 



6 AN OBSERVATIONAL PROBE: LYMAN-SERIES 
ABSORPTION 

One consequence of large ionized bubbles is to reduce or eliminate 
the Ly a "damping wing" absorption by the IGM iFurlanetto et alJ 
2004c). If bubbles remained small throughout reionization, the Lya 
forest would allow no transmission because the damping wings of 
the remaining neutral patches would overlap, even though HII re- 
gions are still highly ionized (Barkana 2002). With large bubbles, 
however, the neutral gas can lie sufficiently far away to render the 
damping wing negligible, opening a new window for studying the 
high-redshift universe: quasar spectra should show interesting fea- 
tures e ven when Xi ~ 0.1-0.25 (see, e.g., Fig. 9 o flFurlanetto et alJ 
l2004d) . 

Our model shows that recombinations strongly affect these 
late stages, so the transmission in this regime constrains how the 
bubbles saturate once the IGM density structure begins to domi- 
nate. To illustrate the power of this approach, we now extend the 
treatment of Furlanetto et al. 1 2004c) to include the recombination 
limit. We will ignore damping wing absorption in the following, 

4 Note that we are being somewhat cavalier with the definition of "Si" 
here; in the FZH04 case it is the bubble filling factor. The true ionized frac- 
tion (by volume) is slightly smaller because of the dense neutral clumps that 
fill a small volume of each bubble. 
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because it is sub-dominant during the epoch of interest. Given a 
model for the bubble topology and iV(A), we can compute the 
probability that any patch of the IGM has an optical depth for tran- 
sition i smaller than n by: 



P{< Ti) = I dm b n b (m b ) ^ 



dAPv(A), 



(21) 



where A max is the maximum density for which r < n. We assume 
ionization equilibrium within each bubble, so the neutral fraction is 

Xenna{T) 



r 



(22) 



where Xe is a correction for (singly-ionized) helium. In the Ap- 
pendix, we show that this assumption breaks down for gas with 
A > Ai, so we always take such regions to be fully neutral. We 
assume an equation of state T = To A 7 , with To = 10 4 K, and we 
approximate a oc T -0 ' 7 . The ionizing rate per hydrogen atom is 



r = / Av e v \i a v 



(23) 



Here the ionization cross-section is approximately a v — 
ctoO/^o) -3 , with (Jo = 6.3 x 1CT 18 cm 2 and hv = 13.6 eV. We 
will write the photon emissivity (by number) per unit frequency as 
(■v = (eo/Vo) (v/i'o)~ rl ~ 1 ', tne tota l emissivity is then er = eo/v- 
We take rj — 3/2 to approximate a starburst spectrum. Finally, be- 
fore reionization, the mean free path is Xi = min(7?b, J? max ); it 
is frequency-independent in the limit that the bubbles are separated 
by walls of neutral hydrogen. (We take a single i? max to repre- 
sent the recombination limit at any instant, which appears to be a 
good assumption unless recombinations are quite rapid.) Thus (c.f., 
Miral da-Escudel2003l) 

XiaoerV 



r 



(24) 



3 + ry 

We can relate the emissivity to the efficiency parameter of FZH04 
by setting er = C,n b xZ (d/coii, g /dt). This assumes e; to be uni- 
form throughout the universe; actually it is slightly larger around 
clusters of sources, but the fluctuations are small at the end of reion- 
ization. We also note that C, is here the instantaneous emissivity and 
should not include any factor accounting for past recombinations. 
Finally, we take the local density field Py(A) to be uncorrelated 
with the bubble topology; when Xi ~ 1, this is probably not a 
bad approximation, though it must eventually be compared to high- 
resolution simulations. 

With these assumptions, we can write the local density as a 
function of observed optical depth, 

qa(10 4 K)] , , , ( 



A(n) 



170 



x C 



V 

3 + V 
d/coii 



h(z) 



Xi 
^Mpc 



dz 



TGP, 



1/(2-0.77) 



(25) 



where tgp,i is the iGunn & Petersonl 1 19651) optical depth for a 
fully neutral medium in transition i; substituting for A max in equa- 
tion 121 i . we can then find the probability that any region of the 
IGM will have a small optical depth in our inhomogeneous bubble 
model. 

Figure [121 shows the resulting probability distributions for the 
Lya, Ly/3, and hy-y transitions. In each panel, the dotted, dot- 
dashed, short-dashed, long-dashed, and solid curves take 7? max = 
10, 20, 30, 60, and 600 Mpc (note that the last is essentially equiv- 
alent to setting -R max = oo). We set Xi = 0.95 at z — 6. The dis- 
tributions peak at a characteristic optical depth determined by the 
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Figure 12. Probability distribution of optical depth for the three Lyman 
transitions, if the universe has Xi = 0.95 at z = 6. The different curves 
assume R max = 10, 20, 30, 60, and 600 Mpc (from left to right). 



convolution of the bubble size distribution with the IGM density 
field. Note how the peak moves relatively little for large 7? max , be- 
cause it can never exceed the characteristic bubble size of FZH04. 
While strong Lya transmission is always relatively rare, the Fig- 
ure shows a sensitive dependence on J? m a X : strong transmission 
relies on deep voids, which are rare in the MHR00 model. Thus 
moving the peak of p(r) effectively moves the cutoff in the opti- 
cal depth distribution and sharply reduces the leftover flux. Similar, 
though somewhat smaller, differences occur for the other Lyman- 
series transitions. The effects are less dramatic here because the 
higher trans itions sample more re presentative parts of the density 
distribution lOh & Furlanettol2005l) . 

Figure fT3l shows how the absorption varies with Xi. We as- 
sume z = 6; the thick curves take i? max = 30 Mpc while the 
thin curves take i? max = 600 Mpc. The dotted, dashed, and solid 
curves have Xi = 0.75, 0.85, and 0.95. In the early stages, 7? max 
makes little difference because i? c har <§C fima X ; hence, the dot- 
ted curves are coincident in each panel. However, as the bub- 
bles grow beyond ii max , more and more of the universe becomes 
recombination-limited and thick and thin curves diverge (see also 
Fig. HOt . At first only the tails of the distributions differ because 
i? max > Rcha.i at Xi — 0.85, but by the end of reionization even 
the peaks separate. Again, the contrast is largest for Lya because 
that samples the most distant tail of the IGM density distribution. 

Thus, the recombination limit ultimately determines the ex- 
tent to which we can use quasar (or gamma-ray burst) spectra to 
study the late stages of reionization. If _R max is small (less than 10 
Mpc, for example), the pre-overlap era will have little or no trans- 
mission because the local ionizing background is small everywhere 
even if the damping wing is suppressed. Hence the Gunn-Peterson 
trough would be complete. If, on the other han d, i? max is comp ara- 
ble to the mean free path at z < 6 (~ 20 Mpc: lFan et all2002f) . we 
expect non-negligible transmission as far back as Xi ~ 0.8. Our 
model therefore makes several testable predictions for the influ- 
ence of recombinations on observations. The key is that the scatter 
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Figure 13. Probability distribution of optical depth for the three Lyman 
transitions. Within each panel, the thick curves assume iJ max = 30 Mpc 
and the thin curves assume i? ma x = 600 Mpc. The dotted, dashed, and 
solid curves are for X{ = 0.75, 0.85, and 0.95, respectively (all at z = 6). 

in apparent overlap increases with J? m ax while the allowed trans- 
mission over that range decreases. The mean free path after overlap 
\(z ~ 6) ~ 20 Mpc likely sets an upper bound on i? max , because 
it normally increases monotonically with time. On the other hand, 
the significant Lyman serie s transmiss ion along the line of sight to 
a z « 6.4 quasar (White et alj|2003|). coupled with the saturated 
absorption to a z ~ 6.3 quasar (Bec ker et alJl200ll) . sets a lower 
bound on i? max . The observed scatter in the "overlap redshift" Az 
should increase as i? max decreases. Figure [Tol shows that taking 
Rm&x{z fa 6-6.4) ~ X(z ~ 6) ~ 20 Mpc is consistent with 
the expected scatter in our model. Moreover, 7? max will not evolve 
strongly over that time (see Fig.Q?), barring strong and sudden evo- 
lution in substructure - e.g., photoevaporation of minihalos. These 
arguments can be tightened both by observations of more high- 
redshift lines of sight and by more detailed simulations of the Lya 
forest at the end of reionization. In general, the rate at which the 
transmission evolves will tell us how the universe passes from the 
"bubble-dominated" topology characteristic of reionization to the 
smoother "web-dominated" topology characteristic of the lower- 
redshift universe. While 21 cm observations will teach us a great 
deal about the early and middle stages of reionizati on (when the 
21 cm signal has a high contrast; FZH04. iFurlanetto et aHz OfMbl 
quasar spectra will likely remain our best probe of the tail end of 
reionization and this crucial overlap epoch. 



7 DISCUSSION 

In this paper, we have studied how HII regions evolve through- 
out reionization. Building on the FZH04 model, in which ion- 
ized bubbles surround clusters of luminous sources, we first de- 
scribed how these bubbles combine and merge in the early and 
middle stages of reionization. We derived bubble progenitor distri- 
butions and (approximate) merger rates analogous to the "extended 
Press-Schechter" quantities commonly used for dark matter halos. 



The distributions emphasize the "inside-out" nature of reionization. 
Large bubbles surrounding overdense regions begin ionizing earlier 
than average, while small bubbles (characteristic of voids) ionize 
later on. They also show that HII regions grow primarily through 
mergers with bubbles slightly larger than the current characteristic 
size. As a consequence, on a local level reionization is more analo- 
gous to "punctuated equilibrium" than "gradualist" evolution, even 
if Xi increases smoothly. Figure|6|shows several examples in which 
bubbles sit quietly for long periods between rapid growth spurts at 
major mergers. 

This behavior changes as we approach the end of reioniza- 
tion. As the bubbles expand, internal recombinations consume an 
ever growing fraction of ionizing photons. In other words, once 
a bubble's size becomes comparable to the mean free path of an 
ionizing photon, it stops growing and becomes a true "Stromgren 
sphere." We showed how this saturation point can be incorporated 
into the FZH04 model by building on the picture presented by 
MHR00. Note that once a significant fraction of the universe sits 
within such bubbles, the bubble sizes do not meaningfully reflect 
variations in the mean free path of ionizing photons. In this regime, 
our model describes large-scale fluctuations in the mean free path. 
This makes interpreting observations somewhat subtle: direct mea- 
sures of the ionized v olumes, such as 21 cm studies of the high- 
redshift universe (e.g. iMadau et alJll997tlCiardi & Madaull2003t 
Zalda rriaga et all2004tlFurlanetto et all2004bl) . will see large con- 
tiguous ionized regions but can tell us little about how the mean free 
path evolves and hence how reionization actually ends. Instead, ob- 
servations of Lyman- series absorption toward high-redshift quasars 
and of Lya-emitting galaxies may remain our best probes of the 
overlap era, because those techniques are much more sensitive to 
the strength of the ionizing background and hence to the mean free 
path of ionizing photons. 

In IGM models similar to MHR00, bubbles saturate only on 
scales > 20 comoving Mpc. However, applying the MHR00 model 
at z > 6 requires extrapolation from their simulation fits. In par- 
ticular, the simulations assumed that the IGM is smooth on the 
Jeans scale of photoionized gas (T ~ 10 4 K). This obviously 
need not be the case for neutral gas at high redshifts. In princi- 
ple the gas may be colder than the CMB and a large fraction could 
condense into minihalos, which act as sinks of ionizing photons. 
Minihalos will significantly reduce the mean free path and hence 
freeze the bubble growth earlier. We considered an extremely sim- 
ple model to illustrate their possible effects; a more detailed look 
must consider the photoevaporation of these ob jects as they are ion- 
ized iBarkana & Loebll999tlshapiro et all2 004). which makes the 
problem time-dependent. Nevertheless, our picture fits with the cur- 
rent observational constraints. iFan et"al]l2002l) measured the mean 
free path of ionizing photons at z ~ 5.5 to be ~ 20 Mpc (using the 
MHR00 model). If the MHR00 model is accurate at high redshifts, 
then the bubbles grow smoothly until this limit is reached when 
Xi > 0.9 (see Fig.^Oj. If the IGM is significantly more clumpy, the 
universe may pass through a long highly-ionized phase in which 
the mean free path grows slowly as, e.g., minihalos photoevapo- 
rate. The true picture awaits a better understanding of the role of 
feedback and IGM heating in suppressing small-scale structure in 
the IGM lOh & HaimaJ2003l) . 

Another effect of our description of inhomogeneous reioniza- 
tion is to increase the recombination rate of the universe above that 
predicted by MHR00 (by up to a factor ~ 4), especially when x% 
is relatively small. This occurs because ionized bubbles originate 
around clusters of galaxies (which in turn occur in overdense re- 
gions) and because, within each bubble, the minimum density to 
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which the IGM must be ionized is larger than for a uniformly il- 
luminated universe. Thus we have found that recombinations ex- 
act a "progressive tax" on the ionizing sources that becomes larger 
around the largest proto-clusters of sources: any such group that 
can blow a large bubble into the IGM will soon have its ionizing 
photons consumed by recombinations, although its more isolated 
neighbors remain effective ionizing agents. This helps to explain 
the simulations of Ciardi et al. 1 2003), who argued that the densest 
regions are the first to begin reionization but the last to "complete" 
it. 

Most importantly, our model shows that the radiation back- 
ground at any given point in the IGM will evolve through a se- 
ries of discrete jumps until it saturates at the mean free path. We 
must therefore be careful in interpreting jumps in the ionizing 
background (or in the IGM transmission) - every point will ex- 
perience a sequence of jumps throughout reionization, with the 
final (and likely largest) jump occurring when a region becomes 
recombination-limited. After this point the ionizing background 
evolves only slowly (see Fig.|6j. Although most regions reach this 
stage over a relatively short time interval (see Fig. llOI . some bub- 
bles grow rapidly and saturate at much earlier times, particularly 
if the limit appears when i?< 10 comoving Mpc. Once a bubble 
reaches this point, reionization is complete for all intents and pur- 
poses within that volume, regardless of the true global Xi (c.f, 
iGnedinlfeood) . and it is obvious that we can no longer treat the 
"overlap epoch" as a single global event. Instead these large bub- 
bles enter a "post-overlap" phase when Si < 1! As a result, we 
expect substantial cosmic variance in the apparent onset of reion- 
ization along different lines of sight. Note that this is a qualitatively 
different source of variance than those described by Wvithe & Loeb 
l2004al) . who attributed differences in the onset of reionization to a 
combination of the light-travel time across bubbles and modulation 
by the underlying density field (similar to the pure FZH04 model). 
While the latter is crucial to the large-scale topology of reioniza- 
tion, the end stages of reionization are actually driven by fluctua- 
tions in the mean free path within large contiguous HII regions and 
not simply by source clustering. In comparison to their model, we 
predict larger scatter in the "overlap" time but smaller scatter in 
the mean free path during that interval, because every region's size 
must be near the Stromgren limit. 

This limit i? max determines the ionizing background within 
the bubbles (and hence the transmission, given a density field and 
equation of state). Thus observations of Lyman-series transmission 
in quasar (or gamma-ray burst) spectra can constrain how substruc- 
ture in the IGM assumes control of the metagalactic radiation field 
and confines the bubbles. If, for example, minihalos limited the 
mean free path to Mpc scales, significant transmission cannot occur 
until reionization is truly complete; on the other hand, density mod- 
els similar to MHR00 predict measurable transmission even when 
i, ~ 0.75. 

Our model also provides a clear and intuitive interpretation 
of the differenc es between the two best-studied lines of sight to 
z > 6 quasars feecker et al]|200ll IWhite et alll2003l) . which go 
much deeper than differences in the endpoint of the Gunn-Peterson 
trough. One shows apparent complete absorption, while the other 
remains highly ionized and has a number of transmission features 
at z > 6. These differences are intriguing because they appear to be 
coherent over ~ 100 comoving Mpc, on which scales density fluc- 
tuations should be miniscule at such high redshifts. In our model, 
the difference arises because overlap occurs somewhat earlier along 
one line of sight. The FZH04 model shows that reionization effec- 
tively amplifies small density fluctuations by distributing the IGM 



gas through a range of bubble sizes. Here we have shown that, dur- 
ing the end stages of reionization, substructure in the IGM further 
modulates the mean free path and hence the expected transmission. 
Thus, at least qualitatively, it appears reasonable that the line of 
sight to the z = 6.28 quasar may pass through a slightly under- 
dense region along which the HII regions have not yet reached the 
Stromgren radius, suppressing the transmission relative to the other 
quasar. In any case, we have shown that such large-scale fluctua- 
tions are a generic feature of reionization and can at least in princi- 
ple be used to infer some key properties of the IGM. 

All that said, our model is approximate and must be im- 
proved through comparisons to observations. FZH04 noted sev- 
eral simplifying assumptions; these included the inhomogeneous 
distribution of sources inside bubbles, the possibility that £ varies 
across galaxies, and the Lagrangian treatment of purely linear fluc- 
tuations. However, the most important were recombinations and 
purely spherical bubbles. We have shown that recombinations do 
affect the bubble growth, but only in the late stages of reioniza- 
tion, helping to validate the FZH04 model. However, all of our 
conclusions are subject to caveats about the unknown density dis- 
tribution within the IGM. That is something that clearly must be 
simulated more carefully in the future if we are to understand the 
full effects of recombinations. It promises to be difficult - follow- 
ing the growth of small-scale structure, including the feedback of 
reionization on it (including the photoevaporation of minihalos; 
Shapiro et al. 2004), while also sampling a representative volume 
of the universe ( > 100 Mpc) - will tax even the promise of Moore's 
law. A hybrid strategy of bootstrapping small boxes, which resolve 
the IGM structure, into larger, lower-resolution boxes may be the 
most viable. 

The other crucial assumption is that of spherical bubbles. We 
argued in Sj2|that this will probably break down when the bub- 
bles fill a large fraction of the universe, because then ionized gas 
can surround isolated islands of neutral gas. In this case the bub- 
ble size may describe the typical separation of such islands. This 
too must be tested in simulations, or at least in nu merical realiza- 
tions of our barrier formalism (e.g. JZahn et all2005l) . However, our 
recombination-based model is less subject to such problems. We 
have seen that, once the bubbles become large, they fragment into 
Stromgren spheres; our model describes the distribution of mean 
free paths. This quantity is not so sensitive to the underlying topol- 
ogy of ionized and neutral gas, because it supposes only the exis- 
tence of large continuous regions of ionized material. 

We thank M. McQuinn for useful discussions. We are grateful 
to the KITP for hospitality during the beginning of the project, un- 
der NSF grant PHY99-07949. SPO also acknowledges NSF grant 
AST-0407084 for support. 
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APPENDIX A: 
MODEL 



A CLOSER LOOK AT THE MHR00 



The interior of a Stromgren sphere during reionization is similar 
to the universe at a lower redshift, when the Lya forest is in ion- 
ization equilibrium. Our physical picture of this situation is heav- 



ily influenced by numerical simulations, which can reproduce the 
z ~ 3 observations remarkably well by assuming a uniform radia- 
tion field r. In this regime, xm oc A. At face value, this is rather 
different from the MHR00 assumption that xui ~ for A < Aj 
and xm ~ 1 for A > A^. Here, we explore whether these two 
views are compatible. In particular, we compute the column density 
distribution function f(Nm) and compare the resulting mean free 
path with the MHR00 ansatz of equation 1 1 21 . We also check the 
validity of the MHR00 assumption that absorption is dominated by 
Lyman limit systems rather than the cumulative photoelectric opac- 
ity of lower column density systems. As a byproduct, we derive a 
number of useful scaling relations. Throughout this Appendix we 
only consider recombination-limited bubbles, where e ~ A. 

It is useful to develop our physical intuition by making some 
(surprisingly accurate) analytic approximations to the MHR00 
model and examining the salient scalings. Let us consider spher- 
ical absorbers with density profiles p oc r~ a . Suppose a radiation 
field r is incident upon the absorbers. We define a Stromgren sur- 
face via Yr 2 oc p 2 r 3 , or r oc r 1/(1 ~ 2a) . Since r oc A _1/q , this 
implies a relation between the radiation field and the overdensity 
up to which the absorber is ionized, T oc A' 20-1 '/ . The mean 
free path of ionizing photons is inversely proportional to the cross- 
section of the absorber, A oc 1/r 2 oc A 2 ^ a . Finally, the recombi- 
nation rate A oc p 2 r 3 oc r ~ 2a+3 oc A (2q ~ 3)/q . For an isothermal 
sphere, a = 2, these relations imply that Y oc A 3 / 2 , A oc A, and 
A oc A 1/2 . Note that self-consistently, Y oc eA w A\ oc A 3/2 . 
By design, the MHR00 density distribution function, equation 
asymptotes to a power law density profile A oc r ~ 3 ^ ,3 ~ 1 ' > 
for A » 1. At z = 3, /3 = 2.35 (derived from a fit to numeri- 



cal simulations), corresponding to p oc 



J , while for z > 6, 



MHR00 assume isothermal density profiles par", implying 
(5 — 2.5. How quickly does the MHR00 density distribution reach 
this limiting behaviour? In their treatment, the radius of an absorber 
isr oc [1 — _FV(Ai)] 1//3 . Hence, we can compute the effective slope 
of the density power law numerically via 



d(logA) 
d(logr) 



d(log(l-^) 1/3 ) 



d(logA) 



(Al) 



which we show in the upper left panel of Figure IA1I for z = 
(3, 6, 9). The density distribution quickly approaches the limiting 
power law density profiles for the overdensities of interest. Thus the 
approximations for isothermal density profiles should be accurate 
for the regime of interest. Because the number density of absorbers 
evolves with redshift (depending on the details of structure forma- 
tion), the normalization has to be obtained numerically. We find for 
z = (3,6,9): 



Amfp 

A 
r_i2 



(41,33,33) (^1) Mpc 



(6.8,4.4,3.0) 



V100V 



= (0.42,2.7,9.1) 



ViooJ 



(A2) 



where A m f p is in physical (proper) units, A is in units of recombi- 
nations per baryon per Hubble time, F_i2 = F/(10 _1 s~ ), and 
we have used equation 1231 to evaluate Y. These power law scal- 
ing relations are compared to the exact expressions (solid lines) in 
Figure lATl at z — 6. They are remarkably accurate and useful for 
quick estimates. 

The MHR00 model postulates that x m w for A < Aj and 
xhi ~ 1 for A > Ai ; this is rather different from the usual ion- 
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Figure Al. Top left: The logarithmic slope a of the density profile, from equation IA1I . The density profiles quickly reach their asymptotic power law scalings. 
The other panels show the proper mean free path A m f p , recombination rate A, and ionization rate r_ 12, as computed numerically (solid lines) and via the 
power law scaling relations (equation lA2l , at z = 6. Comparison at other redshifts show similarly good agreement. 



ization equilibrium assumption (zhi oc A) in Lya forest studies. 
MHR00 argued that their model could accommodate such a scaling 
if one defined Aj such that XHi(Ai) = 0.5. Then, assuming ion- 
ization equilibrium, one obtains virtually the same recombination 
rate A(Ai) as with the original definition of Aj (see their Figure 
2). However, it predicts the wrong scaling between A, and the radi- 
ation field: if aiHl(Ai) =const, then Aj oc V. By contrast, equation 
<A2t implies that Aj oc T 2 ^ 3 . Indeed, from the scaling relations in 
equation <A2t . we have 



apparently contradicting the MHR00 assumption that a;m(Ai) ~ 
O(l). What has gone wrong? 

The key is that we must take the finite absorber sizes - and in 
particular self-shielding - into account. We therefore need t o asso- 
ciate a column density TVhi with a given overdensity A, . ISchavel 
(2001) has shown that many observed properties of the Lya forest 
can be easily understood if the typical size of an overdensity A is 
the local Jeans length Lj. Then, assuming iVm ~ ihiAShLj 
and photo-ionization equilibrium, we obtain 



lEl(Ai 



X e n H a(r) 

r(A0 Al 

(1.4,1.2,1.0) x 10" 



(A3) 



Viooy 



-1/2 



N m = 3.3xl0 17 cm" 2 



A Lyman-limit system which can self-shield from the surrounding 
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radiation field has TVhi ~ 1/co = 1.6xl0 17 cm 2 . From equation 



<A4> . an object enters this regime at 



Alls = (330,61,21) V 



2/3 



(A5) 



for 2 = (3, 6, 9). Note that this expression now has the correct 
power-law scaling. On the other hand, from equation IA2L we have 

At = (180, 51, 23) r 2/ f 2 . (A6) 

Thus, At ~ Alls - From equations <A4> . <A5> . and <A6t . we con- 
clude that the MHR00 model is indeed fully self-consistent, in the 
sense that xhi ~ for A < Aj, and xm ~ 1 for A > Ai. 
This happens because MHR00 calibrated their model to numeri- 
cal simulations, which obviously do contain the relevant physical 
scales for the absorbers. There is an abrupt transition at A, due to 
self-shielding, so that V — *• for A > Ai. Below the threshold, 
ionization equilibrium with T w const and xm oc A are accurate. 

Equation <A4> also allows us to calculate more detailed statis- 
tics of the Lya forest and to compare with observations at z ~ 3. 
In particular, let us calculate the canonical distribution function, 

d 2 A H(z) 1 



HN m ,z) 



dNmdz Ho {1 + zY 



(A7) 



Observe that: 



a, 



Q. g I dAP v (A,z)A 



8nGmu 



N HI Xnlf{Nm,z)dN B 



3H c(l - Y) 
Hence, differentiating both sides, we obtain 



(A8) 



= fit 



-x&N m . 



(A9) 



/(A H i, Z, , . . — ■ . ■ — ■ - ■ v .-, 

dNui 87rGmH 

Using equation <A4t and our power-law scalings for Py(A), we 
have, at z — (3, 6, 9), 

r (-l. 56,6/3,5/3) 



f(N B 



(4.5,4.2,4.7) x 10" 



N 



x r 



-(0.56,2/3,2/3) 



(A10) 



where Nir = Ahi/(10 cm ). The different power law expo- 
nents arise because we have assumed f3 = 2.35 at z = 3 but /3 = 
2.5 at z > 6. In general, f(N m ) oc A^r oc N^ P,3 Y- 20 I 3+1 . 

To test our model, we compare to observations of the 
Lya forest at z ~ 3. The unknown radiation field T-i2 
can be fixed through matchin g the me an flux decrement in 
quasar spectra. Tschirber & Bullockl 1200 3l) find that r_i2(z) = 



10 



-0.24z+0.4 



is a good approximation to the measurements of 
iMcDonald & Miralda-Escudl 1200 ll) in the redshift range z = 
2.3 — 5. If we insert this radiation field into equation <A9t . we 
find excellent agreement with observations. For instance. iHu et alj 



(1995) find a best fit power law f(Nm) = 5.5 x 10 -12 A-^ 1 ' 46 over 
the range N m = 2 x 10 12 - 3 x 10 14 cm" 2 at z = 2.8, whereas 
equation lA9l predicts f(Nm ) = 7.6 x lO" 12 ^ 1 ' 56 . Likewise, 
IStorrie-Lombardi et alJ 1 19941) find an abundance of Lyman-limit 
systems dN^/dz = 3.3 ± 0.6 at z = 4. By comparison, if we 
integrate equation <A9> over all column densities iVm > A^lls = 
1.6 x 10 17 cm -2 , we find dALLs/dz = 3.8. Thus our model seems 
physically reasonable, although the observations suggest rescaling 
the normalization of equation <A9> by a factor f s ~ 0.7 (proba- 
bly due to some ambiguity in the definition of Lj). As a further 
consistency check, equation <A9> predicts a mean free path 



^mfp 



cH- 1 

0- + *) 



dN 
~dz~ 



(76, 20, 6.4) T 2 ^ (A) Mpc.(Al 



1) 



at z = (3,6,9). By comparison, equation <A2t yields A m f p = 
(73, 17, 7.6) T^Mpc, which is in good agreement. 

The sole remaining point to check is the relation between the 
mean free path A m f p of ionizing photons and the distance between 
Lyman-limit systems Alls ; we have assumed that they are roughly 
comparable, A m f p ~ Alls - This amounts to equating Lyman-limit 
systems with contours of overdensity Aj and neglecting the cu- 
mulative photoelectric opacity of lower column-density systems. 
In general, if /(Am) oc A m *, then the effective optical depth 
T c fr oc A m ' s+2 ; as long as 8 < 2, the highest column density 
systems dominate the effective optical depth and hence determine 
the mean free path. More precisely, the absorption probability per 
unit len gth, , for a p hoton at the hydrogen ionization absorption 
edge is I Miralda-Escude 200^): 



1 

An 



r 04 
Jo 



dTT-°(l-e- 



1 



d-n—' 5 A L 



(2.0,2.6) 

Alls 



(A 12) 



for S — 2/3/3 = (1.56, 5/3) at z = 3 and z > 6 respectively. Here 
we have defined r = A?hi/(L6 x 10 17 cm -2 ). Thus, we see that 
for the assumed model the mean free path is indeed comparable to 
the spacing between Lyman-limit systems. At the Lyman edge, it is 
actually a factor ~ 2 smaller, but higher energy photons will reduce 
some of that difference. 

Note that A m f p /ALLS depends only on the slope of column 
density distribution function <5, which in turn depends on the slope 
a of the density profile around absorbers. If they are shallower, 
photoelectric absorption from lower column-density systems be- 
comes more important. In such a picture, absorption is cumulative 
and continuous rather than like a picket-fence. In particular, for 
a < 1.5, 8 > 2 and the opacity is dominated by low column den- 
sity systems (at 8 — 2 itself, we find A m f p /ALLS ~ 0.1). Note 
also that since A oc A ( - 2q_3 ^° ! , recombinations are dominated 
by low overdensities rather than high overdensities for a < 1.5. 
In this case, the recombination rate converges once emission ex- 
ceeds a certain threshold, allowing the Stromgren surface around 
ionized clumps to suddenly move inward and the mean free path 
to inc rease dramatically (see also the discussion in iMadau et alJ 
1999). Thus, the nature of percolation itself depends quite sensi- 
tively on the unknown slope of the density field in dense clumps. 
One could imagine that gas density profiles might be quite differ- 
ent at high redshift, when potential wells are much shallower and 
Jeans smoothing more important. For instance, if there is a period 
of early reionization that raises the entropy of the IGM, gas den- 
sity profiles around minihal os can be much shallower than the dark 
matter profiles I Oh & Haiman 2003). 
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